Introduction
Wheat (Triticum aestivum L.) is widely planted worldwide including arid
and semi-arid regions and is one of the most important food crops (Cossani and Reynolds 2012; Gao et al. 2019). In China, wheat is the
second most widely grown crop after rice (Li et al. 2017). With an increase in human population, it is
estimated that crop production will be doubled by 2050 to meet the food demand
(Li et al.
2017). Therefore, efforts are still needed to increase wheat
yield.
Xinjiang
Uygur Autonomous Region, located in Northwest China, is an important
wheat-producing region. Fruit tree-based intercropping systems are
widely applied by local farmers because that this cultivation method can
prevent soil erosion, decrease water loss, enhance carbon sequestration,
organic carbon, soil nutrient status, and nutrient cycling (Bergeron et al. 2011; Rivest et al. 2009; 2013). Moreover, this method significantly
increases the land equivalent ratio (Yu et al.
2015).
Photosynthesis
and nitrogen fertilizer are two major factors that affect plant production.
More than 90% of grain yield depends on photosynthesis (Makino 2011). During photosynthesis,
photosynthetic pigments absorb light energy and convert it into chemical energy
in the form of carbohydrates, which are needed for plant growth (Wu et al. 2018). The photosynthetic
rate is mainly determined by three aspects, namely, light intensity,
temperature and carbon dioxide concentration. Among these factors, the light
intensity is of great importance. Generally, higher light intensity corresponds
to higher photosynthetic rate. The competition
for light between wheat plants and intercropping fruit trees, which lead to
shading to wheat plants, might decrease photosynthetic rate and leaf area index
in wheat leaves, thus decrease the grain yield of wheat (Kemp and Whingwiri 1980; Mu et al. 2010). For example, apricot-,
walnut-, and jujube-based intercropping systems lead to a decrease of 31.9,
36.2 and 11.3%, respectively in the photosynthetic rate of wheat, compared to
monocultured wheat (Qiao et al. 2019). Shade-grown plant show
a degree of etiolation, which is as another cause of reduced growth and yield
of crop plants (Batool et al. 2020).
Winter
wheat is the second most planted crop in China after rice. Fruit tree-based intercropping
systems are widely favored in several areas of China. However, the productivity
of intercropping crops in many agroforestry systems can be significantly
decreased because of the shading effect from neighboring tree. Many individual
components and pathways operating during plants adaptation to low light stress
have been identified, however the molecular mechanism(s) on the wheat response
to shading are still not well known. There are opposing findings that
shading can increase the leaf area index and photosynthetic rates and delay
leaf senescence and therefore, increase the wheat grain yield (Li et al. 2010;
Xu et al. 2016).
Whether shading increases or decreases the wheat yield depends on the level of
shading (Wang et al. 2003; Guo et al. 2014; Li et al. 2017). However, the molecular mechanisms in response to
different levels of shading in wheat leaves are not well understood.
In
this study, an RNA sequencing (RNA-seq) analysis on wheat flag leaves at
heading stage under different shading conditions (i.e., no shading, and 25 and 75% shading treatments) was performed.
Through a comparative analysis, we surveyed transcriptome changes in wheat
leaves in response to different levels of shading conditions. As a result, we
find some possible genes, which participated in shading response in wheat
leaves, and revealed the possible regulatory mechanism in wheat leaves in
response to different level of shading on transcriptional level.
Materials
and Methods
Experimental site
The winter wheat (Triticum aestivum L. 'Xindong
20') individuals were planted at zepu base (E 77º16',
N 38º10', altitude 1266 m) in Kashi, Xinjiang, China between 2016 and 2017.
Plant materials and experimental design
We
established three treatment groups: the control (sample name: TDC0), 25%
shading group (sample name: TDT1), and 75% shading group (TDT3). The
individuals in the control were grown under natural light. In 25% shading
group, they were shaded by 10% shading was provided at the elongation stage and
25% at the heading stage. In 75% shading group, they were shaded by 30% shading
was done at the elongation stage and 75% at the heading stage. For shading
treatments, individuals were covered by black nylon mesh with different
luminous transmittance. The shade mesh was kept at about 100 cm above the wheat
canopy and was adjusted depending on the plant growth to ensure enough
ventilation. Each shading
treatment included three repeated plots and each plot covered 8 m2
(4 m × 2 m).
Before planting, 150 kg/hm2 urea and 375 kg/hm2
diammonium phosphate were applied as base fertilizers. Then, 150 kg/hm2
urea and 300 kg/hm2 urea were applied at the returning stage and
jointing stage, respectively. The sowing amount of winter wheat was 270 kg/hm2,
and the row spacing was 20 cm. Drip irrigation was carried out 5 times
during the whole season (i.e., at wintering stage, returning stage,
elongation stage, flowering stage, and grain filling stage). All other
cultivation measures were consistent among groups.
Sampling and measurement
The dry weight of wheat plants
was determined at the flowering period and maturity period. The photosynthetic
rate of flag leaves was determined at flowering and filling stages using Li-6400
portable photosynthetic system (Li-Cor, Lincoln, NE, USA) as described previously
(Li et al.
2017). Briefly, flag leaves with uniform fertility progress and uniform
illumination on healthy plants were selected, and the photosynthetic rate
was determined between 9 am and 12 noon on a sunny
day. The artificial light intensity was set as 1200 µmol·m-2·s-1.
At heading stage, samples were collected after three days of complete shading
treatment. Similar-sized ten pieces of flag leaves were collected from each
plot and were mixed as one sample for RNA extraction.
Construction of complementary DNA (cDNA) library and
RNA-sequencing
Total RNA of the flag leaves
was extracted using TRIzol reagent (ThermoFisher, Waltham, MA, USA). RNA degradation and purity
were checked using 1% agarose gels and spectrophotometer (IMPLEN, CA, USA),
respectively. RNA
concentration and integrity were assessed with commercial kits, respectively.
Then, RNA sequencing libraries were generated using 1 μg
RNA from each sample. Qualified cDNA libraries were submitted for Illumina Hiseq platform to perform RNA sequencing.
Data analysis
Raw reads in FASTQ format were
quality-controlled by in-house Perl scripts to obtain clean reads. During this
step, reads containing adapter, ploy-N (>10%), and with low quality (Qphred ≤ 20 in more than 50% bases) were
removed. Q20, Q30 and guanine-cytosine (GC) content were calculated. The clean
reads were then aligned to the wheat reference genome using HISAT
version 2.1.0 (Kim et al. 2015)
with default parameters. Gene expression was estimated as fragments per
kilobase per million (FPKM) base pairs sequenced, which is the most frequently
used method in estimating gene expression levels (Trapnell et al. 2010), using HTSeq (version 0.6.1) with default parameter of -m union (Anders et al. 2015). Generally,
the gene expression with FPKM greater than 0.1 is regarded as a threshold to
determine whether a gene is expressed or not. Correlations between
RNA-sequencing data of each sample were calculated by the Pearson’s
correlation coefficient in R software (version 3.4.2, Free Software Foundation,
Boston, USA). Qualified RNA-sequencing data of samples in each group should
meet the requirement of r2
greater than 0.8.
Differential
expression analyses were performed using the “DESeq2”
package of R (Love et al.
2014). To reduce the false discovery rate, P-values
were adjusted by Benjamini and Hochberg (1995) approach. The
genes fulfilling the thresholds of adjusted P-value
< 0.05 and |log2 fold change (FC)|≥1 were
regarded as differentially expressed genes (DEGs). Venn diagram was drawn to
show overlapped DEGs between each comparison group, and
hierarchical clustering was performed using the “Heatmap”
package of R to compare expression patterns of the DEGs in each group.
Gene
Ontology (GO) database provides gene descriptions and their products, which is
applicable to various species. It classifies functions along three aspects:
biological process (BP), cellular component (CC), and molecular function (MF) (Ashburner et al. 2000; Carbon et al. 2017). GO classification of DEGs was
performed using GOseq (Young et al. 2010) based on hypergeometric distribution.
The
DEGs were mapped to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database
using KOBAS web
server (version 2.0) (Mao et al. 2005)
to examine the dysregulated pathways of DEGs. Corrected P-value <
0.05 was regarded as a threshold to select the most enriched GO terms and KEGG
pathway terms. Transcription factors (TFs) were identified using iTAK software (Zheng et al. 2016)
with default parameters.
Validating DEGs with quantitative reverse
transcription-polymerase chain reaction (qRT-PCR)
The qRT-PCR
was performed to validate RNA sequencing data of 16 DEGs. Total RNA from wheat
leaves collected from the three groups was extracted using TRIzol
reagent. DNaseI was added to remove genome DNA. The
RNA was reverse transcribed into cDNA with a PrimeScript
RT reagent kit (Takara, Dalian, China, Code No. RR037A) and purified with QIAquick PCR
purification kit (Qiagen). The qRT-PCR was
carried out using a light cycler 480 System (Roche) and SYBR GREEN
kit (TaKaRa). The qRT-PCR
procedure was as follows: 95°C for 15 min, 40 cycles at 95°C for 10 s, 58°C for
20 s, and 72°C for 20 s. The primers are listed in Table 1. Relative expression
of genes was calculated by 2–ΔΔCt
method with tubulin as an internal control (Livak and Schmittgen 2001).
Three biological replicates were analyzed and each reaction was performed in
triplicate.
Results
Comparisons of changes on dry weight and photosynthetic
rate under different levels of shading
Under shading conditions, the development of
wheat plants was inhibited. Plants have longer growth period (2 days longer in
TDT1 and 11 days longer in TDT3, respectively) and decreased leaf size. The dry
weight and photosynthetic rate of flag leaves were also determined in each
group. The dry weights of plants in 25% shading group and 75% shading group were
significantly less than the control group at flowering period (P <
0.05). At maturity period, the dry weight of plants in 75%
shading group was significantly less than the control group (P < 0.05), while there was no
significant (P > 0.05) difference
between 25% shading group and the control (Fig. 1A). Similarly, the
photosynthetic rate also showed a continuously decreasing trend among
groups. That is, the photosynthetic rates in 25 and 75% shading group were
significantly lower than the control group (P < 0.05) at flowering
period. However, the photosynthetic rate at the grain period was significantly (P
< 0.05) decreased only in 75% shading group (Fig. 1B).
RNA sequencing data
A total of 113.07 Gb clean reads were obtained after quality control. Q30 of
all samples was larger than 89.65% and GC contents were between 50.63 and
55.54%, indicating high quality of RNA sequencing. Ratios of reads mapped to
the reference
genome of each sample ranged from 88.34 to 90.6% (Table 2). Pearson correlation
analysis showed that the r2
between pair-wise
samples in each group was all larger than 0.9, indicating high repeatability in
biological replicates and high reliability of RNA-sequencing data (Fig. 2A).
Identification of DEGs
Based
on the criteria of adjusted P-value
< 0.05 and log2FC ≥ 1, a total of 12,597 DEGs, including 4,466
upregulated genes and 8,131 downregulated genes were identified between TDT1
and TDC0. A total of 28,349 DEGs were identified between TDT3 and TDC0,
including 12,371 upregulated genes and 15,978 downregulated genes. A total of
5,532 DEGs were identified between TDT3 and TDT1, including 2,296 upregulated
genes and 3,236 downregulated genes (Fig. 2B–D). The total number of DEGs
between TDT3 and TDC0 was 2.25-fold higher than those identified between
TDT1 and TDC0. Venn diagram between TDT1 versus TDC0 and TDT3 versus TDC0
showed that there were 2,591 and 18,343 unique DEGs found in TDT1 and TDT3,
respectively (Fig. Table 1: The primer
sequence used for qRT-PCR
|
Gene bank ID |
GENE |
Sequence (5′-3′) |
Length(bp) |
|
AA086852 |
W-CYP98A1-F |
CGACTACGGCCCGCACTACATC |
179 |
|
|
W-CYP98A1-R |
GTTCCTCACTACCATTGGCTTTCCTT |
|
|
AA0438620 |
W-PPO-F |
GCTGACCGACCTTGACTCCA |
104 |
|
|
W-PPO-R |
CCAGGAACAGCAGCGTCT |
|
|
AA0153830 |
W-At3g21620-F |
ATGGCTAATCAACTCATCCTC |
190 |
|
|
W-At3g21620-R |
AAGGTTCAAGACAAAGCCCAA |
|
|
AA1246570 |
W-ACO1-F |
ACTTCTTCCTCCACTTGCTCT |
101 |
|
|
W-ACO1-R |
ATTTGCATGTATGGCTGTGAGACC |
|
|
AA0864380 |
W-PKL-F |
TTAAGTACCACACGGGCAGT |
157 |
|
|
W-PKL-R |
CCTGAATCCCTACTCTCGTCA |
|
|
AA006795 |
W-ETR1-F |
AGACCGGGAGGCATGTTAGG |
166 |
|
|
W-ETR1-R |
GTGATTTGGTGGCGCAGAG |
|
|
AA1656990 |
W-WHAB1.6-F |
GCAGCAAGGAGAAGCCAGCATC |
214 |
|
|
W-WHAB1.6-R |
GGGAACCACTCGGCGGACAT |
|
|
AA0115100 |
W-Os10g0493600-F |
GCTGGAACCACTTCTACTGCG |
151 |
|
|
W-Os10g0493600-R |
AAGTTACCCTGGGAATCCCTACTA |
|
|
AA0272780 |
W-6-FEH-F |
ATGGCGTATCGCCGTTGCA |
198 |
|
|
W-6-FEH-R |
GCGTAGCCAAGTCGCCCCTC |
|
|
AA0318750 |
W-CIN2-F |
GCGTACGTGTACCGGAGCAGGG |
153 |
|
|
W-CIN2-R |
GACGACGGCGGACGAGGTGT |
|
|
AA1820630 |
W-1-SST-F |
CAAGGACATGGTGAACTGGCG |
127 |
|
|
W-1-SST-R |
GGTGTAGAGCAGGATGACCTGGC |
|
|
AA0822510 |
W-PSB28-F |
CGCTGCCCATGTCCATGCG |
198 |
|
|
W-PSB28-R |
GACCAAGAATTCCGGCGAACAC |
|
|
AA1249400 |
W-HEMH-F |
TCAACCTCTTCGCTGATCC |
119 |
|
|
W-HEMH-R |
GGCATACCCCTCTTTACTCT |
|
|
AA0275620 |
W-PORA-F |
CTTCAGGGCTAGGTTTAGCAAC |
162 |
|
|
W-PORA-R |
CCAAAGACGCCAAGTCCAA |
|
|
AA0130320 |
W-At5g26707-F |
AAAAGCCAAATAGAACACTAG |
141 |
|
|
W-At5g26707-R |
CTAAATGGGCATTCATCTC |
|
|
AA1616020 |
W-CPX-F |
TTGGACAGCGTGAGGCAGTTT |
182 |
|
|
W-CPX-R |
ATCGGAGTTCTTCAAGTCATCAATCAA |
|
|
AA086438 |
W-B-TUBF |
GGGCAAGATGAGCACCAA |
170 |
|
|
W-B-TUBR |
TCCACGAAGTAGGAGGAGTTCTT |
Table 2: The mapped statistics of all
samples
|
Sample name |
TDC0_1 |
TDC0_2 |
TDC0_3 |
TDT3_1 |
TDT3_2 |
TDT3_3 |
TDT1_1 |
TDT1_2 |
TDT1_3 |
|
Total reads |
72518502 |
99514062 |
72946984 |
76479352 |
88167792 |
92834306 |
73170976 |
81721748 |
96336264 |
|
Total mapped |
66846096 (92.18%) |
91920789 (92.37%) |
66857327 (91.65%) |
70580568 (92.29%) |
81864542 (92.85%) |
87292018 (94.03%) |
67642629 (92.44%) |
76507807 (93.62%) |
88704652 (92.08%) |
|
Multiple mapped |
2125572 (2.93%) |
3411006 (3.43%) |
1983283 (2.72%) |
3020043 (3.95%) |
3121944 (3.54%) |
3181544 (3.43%) |
1829799 (2.5%) |
2824753 (3.46%) |
3741568 (3.88%) |
|
Uniquely mapped |
64720524 (89.25%) |
88509783 (88.94%) |
64874044 (88.93%) |
67560525 (88.34%) |
78742598 (89.31%) |
84110474 (90.6%) |
65812830 (89.94%) |
73683054 (90.16%) |
84963084 (88.19%) |
3A). Another Venn diagram
showed that only 1,170 DEGs were continuously differentially expressed in 25
and 75% shading groups (Fig. 3B).
Hierarchical
cluster analysis of DEGs
The FPKM values of DEGs in
each group were used for hierarchical clustering analysis. Gene expression
patterns were divided into six groups (Fig. 4). Cluster A was enriched with
genes showing upregulation in TDT1 and downregulation in TDT3. Cluster B was
enriched with genes showing continuous upregulation in TDT1 and TDT3. Cluster C
was enriched with genes showing downregulation in both TDT1 and TDT3. Cluster D
and E were enriched with genes showing downregulation or upregulation only in
TDT3. Cluster F was enriched with genes showing downregulation in TDT1 and
upregulation in TDT3.
Functional enrichment of DEGs
The GO enrichment analysis of DEGs between TDT1 and TDC0 as well as
between TDT3 and TDC0 were
%20333-342_files/image002.jpg)
Fig. 1: Changes in the dry weight (A),
photosynthetic rate (B) of flag leaves in each group.*P < 0.05 compared with TDC0; # P < 0.05 compared with TDT1
%20333-342_files/image004.jpg)
Fig. 2: Identification of differentially expressed genes (DEGs). (A) Pearson’s correlation analysis between pair-wise samples in each group. (B) Volcano plot of DEGs between TDT1, TDC0. (C) Volcano plot of DEGs between TDT3, TDC0. (D) Volcano plot of DEGs between TDT3, TDT1
performed.
Based on the criteria of corrected P-value
< 0.05, 91 GO terms were significantly enriched by upregulated genes between
TDT1 and TDC0, including 29 cellular component (CC) terms, 17 molecular
function (MF) terms, and 45 biological process (BP) terms. The GO-BP terms
enriched by upregulated genes included cellular homeostasis (Corrected P = 9.34e-06, n = 50), gene expression
(Corrected P =2.64e-05, n = 644) and
cell redox homeostasis (Corrected P =
3.44e-05, n = 38) (Fig. 5A). Besides, a total of 251 GO terms were
significantly enriched by downregulated genes between TDT1 and TDC0, including
15 CC terms, 96 MF terms, and 140 BP terms. The downregulated genes were mainly
enriched in GO-BP terms related with transmembrane transport (Corrected P = 1.29e-22, n = 452), single-organism
process (corrected P = 1.69e-17, n =
2557), phosphate-containing compound metabolic process (corrected P = 1.50e-16, n = 816), and phosphorus
metabolic process (Corrected P =
1.54e-16, n = 818, Fig. 5B).
%20333-342_files/image006.jpg)
Fig. 3: VENN
diagram displays the intersection of differentially expressed genes in each
group. (A) VENN diagram between TDT1 vs TDC0, TDT3 vs TDC0. (B) VENN diagram between TDT1 vs
TDC0, TDT3 vs TDC0, TDT3 vs
TDT1
%20333-342_files/image008.jpg)
Fig. 4: Hierarchical cluster analysis of differentially
expressed genes (DEGs). (A) Heatmap of DEGs. (B) The
six expression patterns of DEGs
Based
on the criteria of corrected P-value
< 0.05, 325 GO terms were significantly enriched by upregulated genes
between TDT3 and TDC0, including 77 CC terms, 77 MF terms, and 171 BP terms.
The GO-BP terms enriched by upregulated genes included ribonucleoprotein
complex biogenesis (Corrected P =
4.29e-21, n = 133), ribosome biogenesis (Corrected P = 2.33e-19, n = 123), cellular nitrogen compound metabolic
process (Corrected P = 3.03e-18, n =
2576), and nitrogen compound metabolic process (Corrected P = 4.25e-17, n = 2693) (Fig. 5C). Besides, a total of 372 GO
terms were significantly enriched by downregulated genes between TDT3 and TDC0,
including 21 CC terms, 137 MF terms, and 214 BP terms. The downregulated genes
were mainly enriched in GO-BP terms related with protein
phosphorylation, including phosphate-containing compound metabolic process
(Corrected P= 2.23e -83, n= 1,785),
phosphorus metabolic process (Corrected P
= 4.26e-83, n = 1,788), phosphorylation (Corrected P=3.09e-81, n = 1,444), and protein phosphorylation
(Corrected P = 1.18e-75, n = 1,339)
(Fig. 5D).
%20333-342_files/image010.jpg)
Fig. 5: Gene ontology of differentially
expressed genes. (A) Gene ontology of the upregulated
genes between TDT1, TDC0. (B) Gene ontology of the downregulated
genes between TDT1, TDC0. (C) Gene ontology of the upregulated
genes between TDT3, TDC0. (D) Gene ontology of the downregulated
genes between TDT3, TDC0
%20333-342_files/image012.jpg)
Fig. 6: KEGG pathway enrichment analysis of differentially expressed genes. (A) Pathway enrichment analysis of the upregulated genes between TDT1, TDC0. (B) Pathway enrichment analysis of the downregulated genes between TDT1, TDC0. (C) Pathway enrichment analysis of the upregulated genes between TDT3, TDC0. (D) Pathway enrichment analysis of the downregulated genes between TDT3, TDC0
KEGG pathway enrichment of DEGs
Based on the criteria of
corrected P-value < 0.05, seven
KEGG terms were significantly enriched by upregulated genes between TDT1 and
TDC0, including ribosome biogenesis in eukaryotes (corrected P = 1.48e-05, n = 41), RNA polymerase
(corrected P = 0.000291, n = 24), cyanoamino acid metabolism (corrected P = 0.000957, n = 1,444), starch and sucrose metabolism (corrected P = 0.0103, n = 145), porphyrin and
chlorophyll metabolism (corrected P =
0.0455, n = 35), circadian rhythm – plant (corrected P = 0.0463, n = 25), and purine metabolism (corrected P = 0.0463, n = 135, Fig. 6A). The
downregulated genes between TDT1 and TDC0 were significantly enriched in ten
KEGG pathways, including plant hormone signal transduction (corrected P = 7.51e-09, n = 141), biosynthesis of
secondary metabolites (corrected P =
0.000217, n = 376), plant-pathogen interaction (corrected P Table 3: Number of unigenes
in each transcription factor family
|
Family |
Number |
|
MYB |
134 |
|
WRKY |
126 |
|
AP2-EREBP |
120 |
|
bHLH |
112 |
|
NAC |
112 |
|
bZIP |
86 |
|
Orphans |
84 |
|
FAR1 |
83 |
|
HB |
69 |
|
C2H2 |
60 |
|
mTERF |
57 |
|
C3H |
56 |
|
G2-like |
52 |
|
GRAS |
46 |
|
AUX/IAA |
43 |
|
HSF |
41 |
|
AB13VP1 |
40 |
|
SET |
39 |
|
SNF2 |
36 |
|
C2C2-Dof |
35 |
|
PHD |
33 |
|
CCAAT |
31 |
|
C2C2-GATA |
30 |
|
GNAT |
30 |
|
ARF |
29 |
|
MADS |
28 |
|
Tify |
27 |
|
TUB |
25 |
|
TRAF |
23 |
Table 4: Gene expression changes of the seven genes in TDT1
versus TDC0, TDT3 versus TDC0 in RNA-seq results
|
Gene |
log2FoldChange
(TDT1
vs TDC0) |
log2FoldChange
(TDT3vs
TDC0) |
|
CYP98A1 |
-2.8049 |
-3.5722 |
|
At3g21620 |
-2.3869 |
-2.0606 |
|
ACO1 |
-6.3004 |
-5.1099 |
|
ETR1 |
-2.2315 |
-2.0732 |
|
1-SST |
-3.7623 |
-4.4975 |
|
PSB28 |
-2.8358 |
-2.8358 |
|
HEMH |
-0.79888 |
-1.2236 |
|
WHAB1.6 |
-1.101 |
-2.9141 |
|
Os10g0493600 |
1.2544 |
1.0797 |
|
6-FEH |
0.67194 |
2.3548 |
|
CIN2 |
1.4464 |
3.4538 |
|
PORA |
2.1397 |
5.1143 |
|
At5g26707 |
2.1447 |
2.2296 |
|
CPX |
0.83205 |
1.70035 |
|
PPO |
2.0957 |
2.4764 |
|
PKL |
1.1797 |
1.9068 |
= 0.000217, n = 89),
phenylpropanoid biosynthesis (corrected P
= 0.000217, n = 86), phenylalanine metabolism (corrected P = 0.00846, n = 64),
glycerophospholipid metabolism (corrected
P = 0.0211, n = 46), metabolic pathways (corrected P = 0.211, n = 649), glutathione metabolism (corrected P = 0.0246, n = 51), nitrogen
metabolism (corrected P = 0.0363, n =
23), and tryptophan metabolism (corrected
P = 0.0476, n = 19, Fig. 6B).
Nine KEGG terms were significantly enriched by upregulated genes
between TDT3 and TDC0, including RNA transport (corrected P = 5.75e-08, n = 145), spliceosome (corrected P = 3.09e-05, n = 156), purine metabolism (corrected P = 0.00012, n = 113), mRNA
surveillance pathway (corrected P =
0.00012, n = 95), RNA degradation (corrected
P = 0.00012, n = 85), homologous recombination (corrected P = 0.00241, n = 52), pyrimidine
metabolism (corrected P = 0.00481, n
= 87), RNA polymerase (corrected P =
0.00712, n = 43), and mismatch repair (corrected P = 0.00860, n = 41, Fig. 6C). Only one KEGG pathway term was
significantly enriched by downregulated genes between TDT3 and TDC0, which was
biosynthesis of secondary metabolites (corrected P = 0.004674, n = 769; Fig. 6D).
TFs involved in wheat development
The iTAK software was used to analyze TFs. We identified 2,015 unigenes
as TFs. The main TF families identified were MYB (134 unigenes),
WRKY (126 unigenes), AP2-EREBP (120 unigenes), bHLH (112 unigenes), NAC (112 unigenes), bZIP (86 unigenes), and FAR1 (83 unigenes) (Table 3).
The qRT-PCR validation of
genes involved in photosynthesis
Since the RNA-Seq data were
high-throughput and might generate some false-discovery results, we used qRT-PCR to validate some important genes in response to
shading treatment. Sixteen DEGs involved in porphyrin and
chlorophyll metabolism, photosynthesis, and galactose metabolism were
selected to be validated, including eight downregulated genes of CYP98A1,
At3g21620, ACO1, ETR1, 1-SST, PSB28, HEMH, and WHAB1.6; and eight upregulated genes of Os10g0493600,
6-FEH, CIN2, PORA, At5g26707, CPX, PPO, and PKL (Table 4). Five of the
eight downregulated genes (At3g21620, ETR1, 1-SST, PSB28, and HEMH) and
five of the eight upregulated genes (Os10g0493600, 6-FEH, CIN2,
PORA, and CPX) were successfully validated by qRT-PCR,
indicating these genes may play important roles in shading conditions (Fig. 7).
Discussion
In the present study, we
sequenced transcriptomes of leaves under different shading levels using Next
Generation RNA-seq to comprehensively characterize gene changes under different
shading levels. A total of 12,597 DEGs were identified between TDT1 and TDC0,
and a total of 28,349 DEGs were identified between TDT3 and
TDC0. In this study, the hierarchical cluster analysis showed that these DEGs
can be divided into six clusters with different gene expression patterns.
Functional results showed that the upregulated genes between TDT1 and TDC0 were
significantly enriched into cellular homeostasis and cell redox homeostasis, and
those between TDT3 and TDC0 were significantly enriched into cellular nitrogen
compound metabolic process and nitrogen compound metabolic process. From this
study, the downregulated genes between TDT1 and TDC0 as well as those between
TDT3 and TDC0 were both significantly enriched in phosphorus metabolic
processes, such as phosphorus metabolic process, protein phosphorylation,
phosphate-containing compound metabolic process and
phosphorylation.
Cell
homeostasis, especially redox homeostasis, is important for the maintenance of
many cellular processes, including the removal of xenobiotics, protection of
protein thiols, maintenance of oxidation-reduction balance, and response to
reactive oxygen species (ROS) (Ayer et al.
2014). Plant cell redox homeostasis is formed as a result of the balance
between ROS and antioxidant interaction, and it acts as a metabolic interface
for signals derived from metabolism and environment (Foyer and Noctor
2005). We assume that shading stress triggers specific networks of biological
processes involved in redox and ionic homeostasis. According to this study,
tens genes
involved in redox homeostasis were upregulated, including several genes
encoding redoxins (GRXC1, GRXC14, and GLRX2),
glutathione S-transferase U10 (GSTU10), and thioredoxin-like proteins (TRL11,
TRL31, CXXS1, CDSP32, CITRX, YLS8, and WCRKC2). Glutaredoxin is a small redox enzyme, which is an
indispensable element of the plant antioxidant system to
keep proteins in their properly reduced states (Rouhier et al. 2006; Meyer et al. 2012). Glutaredoxins
were reported to be activated under several abiotic stresses, such as oxidative
stress and drought (Kanda et al. 2006;
Ding et al. 2019). Usually, a massive ROS
was accumulated in plants when they suffer from abiotic stresses, and the
activation of glutaredoxins could prevent plants from
damage by scavenging ROS (Hasanuzzaman et al. 2013). GSTU10 is a
member of the Tau family of glutathione S-transferases (GSTs), which are plant-specific and
can bind glutathione conjugated fatty acid derivatives (Dixon and Edwards 2009). In a previous study,
the upregulation of GSTU10
was also found in Arabidopsis with
inducible nicotinamide adenine dinucleotide (NAD) overproduction (Petriacq et al. 2012)
and in Arabidopsis infected with Botrytis cinerea (Windram et al.
2012). Thioredoxin-like proteins are homologous with thioredoxins which was
characterized by the presence of six tetratricopeptide repeats in conserved
positions (Lakhssassi et al.
2012). Several members of thioredoxin-like proteins were reported to be
involved in abiotic stress tolerance, such as TTL1, TTL3, and TTL4 %20333-342_files/image014.png)
Fig. 7: Validation of 16 genes by qRT-PCR. *P < 0.05 compared with TDC0; **P < 0.01 compared with TDC0; #P < 0.05 compared with TDT1; ## P < 0.01 compared with TDT1
(Borsani et al. 2002; Rosado et
al. 2006; Lakhssassi et al. 2012). Since
Arabidopsis TTL1 is functioned in
osmotic stress tolerance, a large-scale upregulation of thioredoxin-like
proteins in this study suggests that these proteins may be important for
appropriate responses to shade stress. Generally, wheat development under 25%
shading is adaptively modulated by a steady-state balance between ROS and
phytohormones.
In
this study, in 75%
shading group, the upregulated genes were significantly enriched in nitrogen
compound metabolic process and the downregulated genes
were significantly enriched in phosphorus metabolic processes. Nitrogen and
phosphorus are essential elements in plant growth and they have many roles in
photosynthetic organisms. Plants well-supplemented with nitrogen can produce
higher levels of photosynthetic pigments, and this allows higher electron
transport rates in chloroplasts (Evans and Poorter 2001). However, under
shading, nitrogen allocation is destined to
greater chlorophyll biosynthesis in light-harvesting systems (Hikosaka and Terashima 1995;
1996; Valladares and Niinemets 2008).
This might lead to structural and physiological changes in plants to
acclimatize low irradiance conditions, such as
architectural modifications both in stems and leaves to increase leaf area (Falcioni et al.
2018). These changes tend to reduce the maximum photosynthetic leaf area
capacity, but increase mass-basis photosynthesis (Pearcy et
al. 2005). Therefore, we speculated that wheat survival
under shading might result in the allocation of nitrogen and phosphorus to
acclimatize to low irradiance conditions.
TFs
regulate gene transcriptions. This study showed that several TFs, including
MYB, WRKY, AP2-EREBP, bHLH, and bZIP
family members, were differentially expressed. bZIP TFs are essential in regulating light
and stress signaling, pathogen defense, flower development and seed maturation
(Jakoby et al. 2002). Chang et al. (2019) identified CaLMF as a significant light signaling component mediating
light-regulated camptothecin biosynthesis under
shading treatment. WRKYs are also widely involved in plant responses to biotic
and abiotic stress (Zhu et al. 2013). Similarly, the MYB
gene family plays an important role in regulatory networks that control
development, metabolism and environmental stresses response (Dubos et al. 2010). Wang et
al. (2015) demonstrated that eight MYB genes were apparently induced, and
three genes were repressed under shading treatment.
Conclusion
This study characterized the
transcriptome profile of wheat in response to different shading levels. A total
of 12,597 and 28,349 DEGs were identified in 25 and
75% shading groups, respectively, compared with the control.
Functional enrichment analyses suggested that the upregulated genes in 25%
shading group were mainly involved in regulating cell homeostasis, especially
in redox homeostasis. The upregulated genes in 75%
shading group were involved in nitrogen compound metabolic process
and the downregulated genes of the two groups were both enriched in phosphorus
metabolic processes. Ten of the 16 DEGs involved in porphyrin and chlorophyll
metabolism, photosynthesis, and galactose metabolism were successfully
validated by qRT-PCR. These data provided detailed
transcriptional changes of wheat in response to different shading levels and
contribute to systemically explained gene regulatory networks mediating the
wheat response to shading.
Acknowledgement
This work was financially
supported by the National Key Research and Development Program of China (No.
2016YFD0300110) and the National Natural Science Foundation of China (No.
31560370).
Author’s contribution
Qi Zhao conceived and designed
the experiments; Xiaorong Li and Jianfeng
Li performed the experiments and wrote the paper; Hongzhi Zhang and Lihong Wang
performed the field experiments; Zhong Wang, Yueqiang Zhang and Jia Shi analyzed the data; Zheru Fan contributed materials for use in the experiments.
References
Anders S, P TPyl, W Huber (2015). HTSeq—a
Python framework to work with high-throughput sequencing data. Bioinformatics 31:166‒169
Ashburner M, CA
Ball, JA Blake, D Botstein, HL Butler, JM Cherry, AP Davis, K Dolinski, SS Dwight, JT Eppig, MA
Harris, DP Hill, L Isseltarver, A Kasarskis,
SE Lewis, JC Matese, JE Richardson, M Ringwald, GM
Rubin, G Sherlock (2000). Gene ontology: Tool for the
unification of biology. Nat Genet
25:25‒29
Ayer A, CW Gourlay, IW Dawes (2014).
Cellular redox homeostasis, reactive oxygen species, replicative ageing in
Saccharomyces cerevisiae. FEMS Yeast Res
14:60‒72
Batool S, A Wahid, SMA
Basra, M Shahbaz (2020). Light augments the action of foliar applied plant
growth regulators: Evidence using etiolated maize (Zea mays) seedlings. Intl J
Agric Biol 23:801‒810
Benjamini Y, Y Hochberg (1995). Controlling the
false discovery rate: A practical, powerful approach to multiple testing. J Royal Stat Soc Ser B (Methodological)
57:289‒300
Bergeron M, S
Lacombe, RL Bradley, J Whalen, A Cogliastro, MF
Jutras, P Arp (2011). Reduced soil
nutrient leaching following the establishment of tree-based intercropping
systems in eastern Canada. Agrofor Syst
83:321‒330
Borsani O, J Cuartero, V Valpuesta, MA Botella (2002). Tomato tos1 mutation identifies a gene
essential for osmotic tolerance, abscisic acid sensitivity. Plant J 32:905‒914
Carbon S, J Chan, R
Kishore, R Lee, HM Muller, D Raciti, KV Auken, PW Sternberg (2017). Expansion of the Gene Ontology
knowledgebase, resources. Nucl Acids Res
45:D331‒D338
Chang CH, ZW Liu, YY Wang, ZH Tang, F Yu (2019). A bZIP transcription factor, CaLMF, mediated light-regulated camptothecin
biosynthesis in Camptotheca acuminata.
Tree Physiol
39:372‒380
Cossani CM, MP Reynolds
(2012). Physiological traits for improving heat tolerance in wheat. Plant
Physiol 160:1710‒1718
Ding SC, FY He, WL
Tang, HW Du, HW Wang (2019). Identification of maize CC-type glutaredoxins that are
associated with response to drought stress. Genes
Basel 10; Article 610
Dubos C, R Stracke, E Grotewold, B Weisshaar, C Martin,
L Lepiniec (2010). MYB transcription factors in
Arabidopsis. Trends Plant Sci 15:573‒581
Evans JR, H Poorter (2001). Photosynthetic
acclimation of plants to growth irradiance: the relative importance of specific
leaf area, nitrogen partitioning in maximizing carbon gain. Plant Cell Environ 24:755‒767
Falcioni R, T Moriwaki,
E Benedito, CM Bonato, LA
de Souza, WC Antunes (2018). Increased
gibberellin levels enhance the performance of light capture efficiency in
tobacco plants, promote dry matter accumulation. Exp Plant Physiol 30:235‒250
Foyer CH, G Noctor (2005). Redox homeostasis, antioxidant signaling: a
metabolic interface between stress perception, physiological responses. Plant Cell 17:1866‒1875
Gao JW, QC Luo, CJ
Sun, H Hu, F Wang, ZW Tian, D Jiang, WW Cao, TB Dai (2019). Low nitrogen
priming enhances photosynthesis adaptation to water-deficit stress in winter
wheat (Triticum aestivum
L.) Seedlings. Front Plant Sci 10; Article 818
Hasanuzzaman M, K Nahar, MM Alam, R Roychowdhury, M Fujita
(2013). Physiological, biochemical, molecular mechanisms of heat stress
tolerance in plants. Intl J Mol Sci 14:9643‒9684
Hikosaka K, I Terashima (1995). A model of the acclimation of photosynthesis in the leaves of C3 plants
to sun, shade with respect to nitrogen use. Plant Cell Environ 18:605‒618
Hikosaka K, I Terashima (1996). Nitrogen partitioning among photosynthetic
components, its consequenses in sun, shade plants. Funct Ecol 10:335‒343
Jakoby M, B Weisshaar,
W Dröge-Laser, J Vicente-Carbajosa,
J Tiedemann, T Kroj, F Parcy
(2002). bZIP transcription factors in Arabidopsis. Trends
Plant Sci 7:106‒111
Kanda M, Y Ihara, H
Murata, Y Urata, T Kono, J Yodoi,
S Seto, K Yano, T Kondo (2006). Glutaredoxin
modulates platelet-derived growth factor-dependent cell signaling by regulating
the redox status of low molecular weight protein-tyrosine phosphatase. J
Biol Chem 281:28518‒28528
Kemp D, E Whingwiri (1980). Effect
of tiller removal, shading on spikelet development, yield components of the
main shoot of wheat, on the sugar concentration of the ear, flag leaf. Aust J Plant Physiol
7:501‒510
Kim D, B Langmead,
SL Salzberg (2015). HISAT: A
fast spliced aligner with low memory requirements. Nat Methods 12:357‒360
Lakhssassi N, VG Doblas, A Rosado, AE Del Valle, D Pose, AJ Jimenez, AG
Castillo, V Valpuesta, O Borsani,
MA Botella (2012). The Arabidopsis tetratricopeptide thioredoxin-like
gene family is required for osmotic stress tolerance, in male sporogenesis. Plant Physiol 158:1252‒1266
Li HW, D Jiang, B Wollenweber, TB Dai, WW Cao (2010). Effects of shading on morphology, physiology, grain yield of winter
wheat. Eur J Agron
33:267‒275
Li Q, SF Zhong, SF Sun, SA Fatima, M Zhang, WQ
Chen, QL Huang, SW Tang, PG Luo (2017). Differential effect of whole-ear
shading after heading on the physiology, biochemistry, yield index of
stay-green, non-stay-green wheat genotypes.
PLoS One 12:Article
e0171589
Livak J, TD Schmittgen (2001). Analysis of relative gene expression
data using real-time quantitative PCR, the 2(-Delta Delta
C(T)) method. Methods 25:402‒408
Love MI, W Huber, S
Anders (2014). Moderated estimation of fold change, dispersion for RNA-seq data
with DESeq2. Genome Biol 15; Article 550
Makino A (2011).
Photosynthesis, grain yield, nitrogen utilization in rice, wheat. Plant
Physiol 155:125‒129
Mao XZ, T Cai, JG Olyarchuk, LP Wei (2005).
Automated genome annotation, pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics 21:3787‒3793
Meyer Y, C Belin, V Delormehinoux, J Reichheld, C Riondet
(2012). Thioredoxin, glutaredoxin systems in plants:
molecular mechanisms, crosstalks, functional significance. Antioxid Redox Signal 17:1124‒1160
Mu H, D Jiang, B Wollenweber, TB Dai, QJing, WW
Cao (2010). Long-term low radiation
decreases leaf photosynthesis, photochemical efficiency, grain yield in winter
wheat. J Agron
Crop Sci 196:38‒47
Pearcy RW, H Muraoka, F Valladares (2005). Crown
architecture in sun, shade environments: assessing function, trade-offs with a
three-dimensional simulation model. New
Phytol 166:791‒800
Petriacq P, LD Bont, J Hager, L Didierlaurent, C
Mauve, F Guerard, G Noctor, S Pelletier, J Renou, G Tcherkez, B Gakiere (2012). Inducible NAD overproduction in Arabidopsis
alters metabolic pools, gene expression correlated with increased salicylate
content, resistance to Pst-AvrRpm1. Plant
J 70:650‒65
Qiao X, LH Sai, X Chen, LH Xue, JJ Lei (2019). Impact of fruit-tree shade intensity on
the growth, yield, quality of intercropped wheat. PLoS One14; Article e0203238
Rivest D, M Lorente, A Olivier, C Messier, (2013). Soil biochemical
properties, microbial resilience in agroforestry systems: Effects on wheat growth
under controlled drought, flooding conditions. Sci Total Environ 463–464:51‒60
Rivest D, A Cogliastro,
A Olivier (2009). Tree-based intercropping systems increase growth, nutrient
status of hybrid poplar: A case study from two Northeastern American
experiments. J Environ Manage 91:432‒440
Rosado A, AL Schapire, RA Bressan, AL Harfouche, PM
Hasegawa, V Valpuesta, MA Botella (2006). The
Arabidopsis tetratricopeptide repeat-containing protein TTL1 is required for
osmotic stress responses, abscisic acid sensitivity. Plant Physiol 142:1113‒1126
Rouhier N, J Couturier, J Jacquot (2006). Genome-wide analysis of plant glutaredoxin systems. J
Exp Bot 57:1685‒1696
Trapnell C, BA Williams, G Pertea, A Mortazavi, G Kwan, MJV Baren,
SL Salzberg, BJ Wold, L Pachter (2010). Transcript assembly, quantification by
RNA-Seq reveals unannotated transcripts, isoform switching during cell
differentiation. Nat Biotechnol
28:511‒515
Valladares F, U Niinemets
(2008). Shade tolerance, a key plant
feature of complex nature, consequences. Annu Rev Ecol Evol Syst
39:237‒257
Wang FQ, YF Suo, H Wei, MJ Li, CX Xie, LN Wang, XJ Chen, ZY Zhang (2015). Identification,
characterization of 40 isolated Rehmannia glutinosa MYB
family genes, their expression profiles in response to shading, continuous
cropping. Intl J Mol Sci 16:15009‒15030
Wang Z, Y Yin, M He, Y
Zhang, S Lu, Q Li, S Shi (2003). Allocation
of photosynthates, grain growth of two wheat cultivars with different potential
grain growth in response to pre-, post-anthesis shading. J Agron Crop Sci
189:280‒285
Windram OP, P Madhou, S Mchattie, C Hill, R
Hickman, EJ Cooke, DJ Jenkins, CA Penfold, L Baxter, E Breeze, SJ Kiddle, J Rhodes, S Atwell, DJ Kliebenstein,
Y Kim, O Stegle, KM Borgwardt,
C Zhang, A Tabrett, R Legaie,
JD Moore, B Finkenstadt, DL Wild, A Mead, DA Rand, J
Beynon, S Ott, V Buchananwollaston, KJ Denby (2012). Arabidopsis defense against Botrytis cinerea: Chronology, regulation deciphered by
high-resolution temporal transcriptomic analysis. Plant Cell 24:3530‒3557
Wu HY, NR Shi, XY An, C Liu, HF Fu, L Cao, Y
Feng, DJ Sun, LL Zhang (2018). Candidate genes for yellow leaf color in common
wheat (Triticum aestivum
L.), major related metabolic pathways according to transcriptome profiling. Intl J Mol Sci 19; Article e1594
Xu CL, HB Tao, P
Wang, ZL Wang (2016). Slight shading
after anthesis increases photosynthetic productivity, grain yield of winter
wheat (Triticum aestivum L.) due to the
delaying of leaf senescence. J Integr Agric 15:63‒75
Young MD, MJ Wakefield, GK Smyth, A Oshlack (2010). Gene ontology
analysis for RNA-seq: accounting for selection bias. Genome Biol 11;
Article R14
Yu Y, TJ Stomph, D Makowski, Wd Werf
(2015). Temporal niche differentiation
increases the land equivalent ratio of annual intercrops: A meta-analysis. Field Crops Res 184:133‒144
Zheng Y, C Jiao, HH Sun, HG Rosli, MA Pombo, PF Zhang, M Banf, XB Dai, GB Martin, JJ Giovannoni,
PX Zhao, SY Rhee, ZJ Fei (2016). iTAK: A program for
genome-wide prediction, classification of plant transcription factors,
transcriptional regulators, protein kinases.
Mol Plant 9:1667‒1670
Zhu XL, SW Liu, C
Meng, LM Qin, LN Kong, GM Xia (2013). WRKY transcription factors in wheat, their induction by biotic, abiotic
stress. Plant Mol Biol Rep
31:1053‒1067